# Compositional time series analysis of state budgets
# using seemingly unrelated regression of additive logratios
# with optional lags and state fixed effects
#
# Helper functions (sourced in)
#
# Christpher Adolph, Christian Breunig, and Chris Koski
# March 2018

# R functions

setx.make <- function(nscen,x,cs) {
  x <- na.omit(x)
  xmean <- apply(x,2,"mean")
  xscen <- list(xhyp=NULL,xhyp0=NULL)
  xscen$xhyp <- matrix(data=xmean,nrow=nscen,ncol=ncol(x),byrow=TRUE)
  xscen$xhyp <- cbind(xscen$xhyp,matrix(data=1/length(cs),nrow=nscen,ncol=length(cs)))
  xscen$xhyp0 <- xscen$xhyp
  colnames(xscen$xhyp) <- c(names(x),cs)
  colnames(xscen$xhyp0) <- c(names(x),cs)
  xscen
}

setx.change <- function(xscen,covname,hyp,hyp0=NULL,scen=1) {
  xscen$xhyp[scen,covname] <- hyp
  if (!is.null(hyp0))
    xscen$xhyp0[scen,covname] <- hyp0
  xscen
}

setx.name <- function(xscen,name,scen=1) {
  xscen$name[scen] <- name
  xscen
}

getvar <- function(getvarstring) {
  returnvalue <- eval(parse(text=getvarstring[1]))
  if (length(getvarstring)>1)
    for (i in 2:length(getvarstring)) {
      returnvalue <- cbind(returnvalue,eval(parse(text=getvarstring[i])))
    }
  returnvalue <- as.matrix(returnvalue)
  colnames(returnvalue) <- getvarstring
  returnvalue <- as.data.frame(returnvalue)
  returnvalue
}
  

lag.panel <- function(x,c,t,lagnum) {
  x <- as.matrix(x)
  listc <- unique(c);
  numtper <- nrow(unique(t));
  outmat <- matrix(nrow=nrow(x),ncol=ncol(x))
  
  for (i in 1:nrow(listc)) {
    xc <- as.matrix(x[c==listc[i,],])
    outmat[((i-1)*numtper+1+lagnum):(i*numtper),] <-
      xc[1:(nrow(xc)-lagnum),]
  }
  outmat
}

postest <- function(xhyp,
                     xhyp0=NULL,
                     b=NULL,
                     vc=NULL,
                     sims=1000,
                     ci=0.95,
                     simbetas=NULL,
                     transform="none", # exp, prop, propexp
                     dofirstdiff=TRUE
                     ) {


xhyp <- as.matrix(xhyp)
if (!is.null(xhyp0))
  xhyp0 <- as.matrix(xhyp0)
if (is.null(xhyp0)) fd <- FALSE else fd <- TRUE

# Create storage matrices
os <- matrix(data=1,nrow=sims,ncol=1)
simy <- matrix(data=NA,nrow=sims,ncol=1)
xhyp <- matrix(data=xhyp,nrow=sims,ncol=nrow(xhyp),byrow=TRUE)
               #t(xhyp)*os
if (fd==TRUE) {
  simy0 <- matrix(data=NA,nrow=sims,ncol=1)
  xhyp0 <- matrix(data=xhyp0,nrow=sims,ncol=nrow(xhyp0),byrow=TRUE)
}

if (is.null(simbetas))
  simbetas <- mvrnorm(n = sims, mu=b, Sigma=vc)

simysave <- simy0save <- NULL

for (i in 1:sims) {
  simy[i] <- simbetas[i,,drop=FALSE]%*%t(xhyp[i,,drop=FALSE])
  if (fd) simy0[i] <- simbetas[i,,drop=FALSE]%*%t(xhyp0[i,,drop=FALSE])
}

# Transformations of response
if (transform=="exp") {
  if (fd)
    if (dofirstdiff) simmu <- exp(simy) - exp(simy0)
    else { simmu0 <- exp(simy0)
           simmu <- exp(simy)
         }
  else simmu <- exp(simy)
}

if (transform=="propexp") {
  if (fd)
    if (dofirstdiff) simmu <- (exp(simy) - exp(simy0))/exp(simy0)
    else { simmu0 <- exp(simy0)
           simmu <- exp(simy)
         }
  else simmu <- exp(simy)
}

if (transform=="prop") {
  if (fd)
    if (dofirstdiff) simmu <- (simy - simy0)/simy0
    else { simmu0 <- simy0
           simmu <- simy
         }
  else simmu <- simy
}

if (transform=="none") {
  if (fd)
    if (dofirstdiff) simmu <- simy - simy0
    else { simmu0 <- simy0
           simmu <- simy
         }
  else simmu <- simy
}
simysave <- cbind(simysave,simmu)
if (fd&&!dofirstdiff)
  simy0save <- cbind(simy0save,simmu0)
  
# Iterate through periods
simysort <- sort(simmu)
yhyp <- mean(simmu)
seyhyp <- sd(simmu)
upyhyp <- simysort[trunc(sims*((1+ci)/2))]
lwyhyp <- simysort[trunc(sims*((1-ci)/2))]

# Construct and save output  
dataout <- list(est=yhyp,
                lower=lwyhyp,
                upper=upyhyp,
                se=seyhyp,
                simbetas=simbetas,
                simy=simysave,
                simy0=simy0save
                )
dataout
}



postestd <- function(xhyp,
                     xhyp0=NULL,
                     b=NULL,
                     vc=NULL,
                     t=1,
                     sims=1000,
                     ci=0.95,
                     lagcols=1,
                     simbetas=NULL,
                     transform="none", # exp, prop, propexp
                     dofirstdiff=TRUE
                     ) {


xhyp <- as.matrix(xhyp)
if (!is.null(xhyp0))
  xhyp0 <- as.matrix(xhyp0)
if (is.null(xhyp0)) fd <- FALSE else fd <- TRUE

# Create storage matrices
yhyp <- matrix(data=0,nrow=t,ncol=1)
seyhyp <- matrix(data=0,nrow=t,ncol=1)
upyhyp <- matrix(data=0,nrow=t,ncol=1)
lwyhyp <- matrix(data=0,nrow=t,ncol=1)
os <- matrix(data=1,nrow=sims,ncol=1)
simy <- matrix(data=NA,nrow=sims,ncol=1)
xhyp <- matrix(data=xhyp,nrow=sims,ncol=nrow(xhyp),byrow=TRUE)
               #t(xhyp)*os
laglist <- matrix(nrow=sims,ncol=length(lagcols))
if (fd==TRUE) {
  simy0 <- matrix(data=NA,nrow=sims,ncol=1)
  xhyp0 <-  matrix(data=xhyp0,nrow=sims,ncol=nrow(xhyp0),byrow=TRUE)
  laglist0 <- matrix(nrow=sims,ncol=length(lagcols))
}

for (i in 1:length(lagcols)) {
  laglist[,i] <- xhyp[,lagcols[i]]
  if (fd) laglist0[,i] <- xhyp0[,lagcols[i]]
}

if (fd) {
  xhyp0.lag <- xhyp0
  xhyp.lag <- xhyp0
} else {
  xhyp.lag <- xhyp
}

if (is.null(simbetas))
  simbetas <- mvrnorm(n = sims, mu=b, Sigma=vc)

simysave <- simy0save <- NULL
countt <- 1;
while (countt<=t) {
  # Set up lags
  for (i in 1:length(lagcols)) {
    xhyp[,lagcols[i]] <- laglist[,i]
    if (fd) xhyp0[,lagcols[i]] <- laglist0[,i]
  }
  for (i in 1:sims) {
    simy[i] <- simbetas[i,,drop=FALSE]%*%t(xhyp[i,,drop=FALSE])
    if (fd) simy0[i] <- simbetas[i,,drop=FALSE]%*%t(xhyp0[i,,drop=FALSE])
  }

  # Transformations of response
  if (transform=="exp") {
    if (fd)
      if (dofirstdiff) simmu <- exp(simy) - exp(simy0)
      else { simmu0 <- exp(simy0)
             simmu <- exp(simy)
           }
    else simmu <- exp(simy)
  }

  if (transform=="propexp") {
    if (fd)
      if (dofirstdiff) simmu <- (exp(simy) - exp(simy0))/exp(simy0)
      else { simmu0 <- exp(simy0)
             simmu <- exp(simy)
           }
    else simmu <- exp(simy)
  }

  if (transform=="prop") {
    if (fd)
      if (dofirstdiff) simmu <- (simy - simy0)/simy0
      else { simmu0 <- simy0
             simmu <- simy
           }
    else simmu <- simy
  }

  if (transform=="none") {
    if (fd)
      if (dofirstdiff) simmu <- simy - simy0
      else { simmu0 <- simy0
             simmu <- simy
           }
    else simmu <- simy
  }
  simysave <- cbind(simysave,simmu)
  if (fd&&!dofirstdiff)
    simy0save <- cbind(simy0save,simmu0)
  
  # Iterate through periods
  simysort <- sort(simmu)
  yhyp[countt] <- mean(simmu)
  seyhyp[countt] <- sd(simmu)
  upyhyp[countt] <- simysort[trunc(sims*((1+ci)/2))]
  lwyhyp[countt] <- simysort[trunc(sims*((1-ci)/2))]
  laglist <- cbind(simy,laglist)[,1:ncol(laglist),drop=FALSE]
  xhyp.lag <- xhyp
  if (fd) {
    xhyp0.lag <- xhyp0
    laglist0 <- cbind(simy0,laglist0)[,1:ncol(laglist0),drop=FALSE]
  }
  countt <- countt+1
}

# Construct and save output  
dataout <- list(est=yhyp,
                lower=lwyhyp,
                upper=upyhyp,
                se=seyhyp,
                simbetas=simbetas,
                simy=simysave,
                simy0=simy0save
                )
dataout
}

# Make pretty ticks & limits
autotic <- function(data,maxtics=12,lower=NA,upper=NA,ntics=6) {
  range <- limits <- matrix(nrow=ncol(data),ncol=2)
  at <- matrix(nrow=ncol(data),ncol=maxtics)
  for (i in 1:ncol(data)) {
    limits[i,1] <- range[i,1] <- min(data[,i])
    factor <- 1
    while (abs(limits[i,1])<1) {
      factor <- factor*10
      limits[i,1] <- limits[i,1]*factor
    }
    limits[i,1] <- floor(range[i,1]*factor)/factor

     
    limits[i,2] <- range[i,2] <- max(data[,i])
    factor <- 1
    while (abs(limits[i,2])<1) {
      factor <- factor*10
      limits[i,2] <- limits[i,2]*factor
    }
    limits[i,2] <- ceiling(range[i,2]*factor)/factor
  }

  if (!is.na(lower))
    limits[,1] <- rep(lower,nrow(limits))

  if (!is.na(upper))
    limits[,2] <- rep(upper,nrow(limits))
  
  at <- matrix(nrow=nrow(limits),ncol=maxtics)
  for (i in 1:nrow(at))
    at[i,1:length(pretty(limits[i,],n=(ntics-1)))] <- pretty(limits[i,],n=(ntics-1))

  for (i in 1:nrow(limits)) {
    limits[i,1] <- min(limits[i,1],na.omit(at[i,]))
    limits[i,2] <- max(limits[i,2],na.omit(at[i,]))
  }
  
  value <- list(limits=limits,
                at=at)
  value
}
  


printnicetable <- function(b,vc,neqs,digits=5) {
  nvars <- length(b)/neqs
  b <- matrix(b,ncol=neqs,nrow=nvars,byrow=F)
  se <- sqrt(diag(vc))
  se <- matrix(se,ncol=neqs,nrow=nvars,byrow=F)
  out <- matrix(ncol=neqs,nrow=2*nvars)
  format(b,digits=digits)
  format(se,digits=digits)
  for (i in 1:neqs) {
    for (j in 1:nvars) {
      if (b[j,i]<0)
        out[((j-1)*2+1),i] <- paste("=\"",substr(as.character(b[j,i]),1,digits+1),"\"",sep="")
      else
        out[((j-1)*2+1),i] <- paste("=\"",substr(as.character(b[j,i]),1,digits),"\"",sep="")
      if (se[j,i]<0)
        out[((j-1)*2+2),i] <- paste("=\"(",substr(as.character(se[j,i]),1,digits+1),")\"",sep="")
      else
        out[((j-1)*2+2),i] <- paste("=\"(",substr(as.character(se[j,i]),1,digits),")\"",sep="")
    }
  }
  out
}



postestdnosim <- function(xhyp,
                          xhyp0=NULL,
                          b=NULL,
                          t=1,
                          lagcols=1,
                          transform="none", # exp, prop, propexp
                          dofirstdiff=TRUE
                          ) {


xhyp <- as.matrix(xhyp)
if (!is.null(xhyp0))
  xhyp0 <- as.matrix(xhyp0)
if (is.null(xhyp0)) fd <- FALSE else fd <- TRUE

n <- nrow(xhyp)

# Create storage matrices
yhyp <- matrix(data=NA,nrow=n,ncol=t)
os <- matrix(data=1,nrow=n,ncol=1)
laglist <- matrix(nrow=n,ncol=length(lagcols))

if (fd==TRUE) {
  yhat0 <- matrix(data=NA,nrow=n,ncol=1)
  laglist0 <- matrix(nrow=n,ncol=length(lagcols))
}

for (i in 1:length(lagcols)) {
  laglist[,i] <- xhyp[,lagcols[i]]
  if (fd) laglist0[,i] <- xhyp0[,lagcols[i]]
}

if (fd) {
  xhyp0.lag <- xhyp0
  xhyp.lag <- xhyp0
} else {
  xhyp.lag <- xhyp
}

yhatsave <- yhat0save <- NULL
countt <- 1;
while (countt<=t) {
  # Set up lags
  for (i in 1:length(lagcols)) {
    xhyp[,lagcols[i]] <- laglist[,i]
    if (fd) xhyp0[,lagcols[i]] <- laglist0[,i]
  }

  yhat <- xhyp%*%b
  if (fd) yhat0 <- xhyp0%*%b


  # Transformations of response
  if (transform=="exp") {
    if (fd)
      if (dofirstdiff) simmu <- exp(yhat) - exp(yhat0)
      else { simmu0 <- exp(yhat0)
             simmu <- exp(yhat)
           }
    else simmu <- exp(yhat)
  }

  if (transform=="propexp") {
    if (fd)
      if (dofirstdiff) simmu <- (exp(yhat) - exp(yhat0))/exp(yhat0)
      else { simmu0 <- exp(yhat0)
             simmu <- exp(yhat)
           }
    else simmu <- exp(yhat)
  }

  if (transform=="prop") {
    if (fd)
      if (dofirstdiff) simmu <- (yhat - yhat0)/yhat0
      else { simmu0 <- yhat0
             simmu <- yhat
           }
    else simmu <- yhat
  }

  if (transform=="none") {
    if (fd)
      if (dofirstdiff) simmu <- yhat - yhat0
      else { simmu0 <- yhat0
             simmu <- yhat
           }
    else simmu <- yhat
  }
  yhatsave <- cbind(yhatsave,simmu)
  if (fd&&!dofirstdiff)
    yhat0save <- cbind(yhat0save,simmu0)
  
  # Iterate through periods
  yhyp[,countt] <- simmu
  laglist <- cbind(yhat,laglist)[,1:ncol(laglist),drop=FALSE]
  xhyp.lag <- xhyp
  if (fd) {
    xhyp0.lag <- xhyp0
    laglist0 <- cbind(yhat0,laglist0)[,1:ncol(laglist0),drop=FALSE]
  }
  countt <- countt+1
}

# Construct and save output  
dataout <- list(est=yhyp,
                yhat=yhatsave,
                yhat0=yhat0save
                )
dataout
}
